#include <Algorithm/projective_recon.h>
#include <Mesh/libcell_iterator.h>
#include <Matrix/LB_Matrix.h>
#include <Math/LB_Quaternions.h>
#include<libcell_macros.h>
#include <Algorithm/threeD_elementary_geometry_tool_library.h>
#include<tool/libcell_tools_algorithm.h>
#define DEFAULT_THRESHOLD  1e-6

#define quote lib_cell_quote

static int is_two_cell_maybe_inter(template_c* c1,template_c* c2)
{
    double min1[3]={c1->vertices[0]->point[0],c1->vertices[0]->point[1],c1->vertices[0]->point[2] },
        max1[3]={c1->vertices[0]->point[0],c1->vertices[0]->point[1],c1->vertices[0]->point[2]};

    for(int i=1;i<c1->vertices_size;i++)
    {
        for(int j=0;j<3;j++)
        {
            c1->vertices[i]->point[j]<min1[j]?(min1[j]=c1->vertices[i]->point[j]):(c1->vertices[i]->point[j]>max1[j]?(max1[j]=c1->vertices[i]->point[j]):i=i);

        }
    }
    double min2[3]={c2->vertices[0]->point[0],c2->vertices[0]->point[1],c2->vertices[0]->point[2] },
        max2[3]={c2->vertices[0]->point[0],c2->vertices[0]->point[1],c2->vertices[0]->point[2]};

    for(int i=1;i<c2->vertices_size;i++)
    {
        for(int j=0;j<3;j++)
        {
            c2->vertices[i]->point[j]<min2[j]?(min2[j]=c2->vertices[i]->point[j]):(c2->vertices[i]->point[j]>max2[j]?(max2[j]=c2->vertices[i]->point[j]):i=i);

        }
    }
    if(min2[0]>max1[0]||min2[1]>max1[1]||min2[2]>max1[2])
    {
        return 0;
    }
    if(max2[0]<min1[0]||max2[1]<min1[1]||max2[2]<min1[2])
    {
        return 0;
    }
    return 1;
}



// static inline int  is_half_sided_co_directional(Mesh*m,template_f*f)
// {
//     if(f->halffaces[0]->vertices[0]==f->halffaces[1]->vertices[1])
//     {
//         return 1;
//     }
//     return 0;
// }
static int is_two_cell_maybe_inter2d(template_c* c1,template_c* c2)
{
    Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(c1->vertices[0]->prop);

    double min1[2]={gpvp->p[0],gpvp->p[1] },
        max1[2]={gpvp->p[0],gpvp->p[1]};


    for(int i=1;i<c1->vertices_size;i++)
    {
        gpvp=(Greedy_Projection_Vert_Prop*)(c1->vertices[i]->prop);
        for(int j=0;j<2;j++)
        {
            gpvp->p[j]<min1[j]?(min1[j]=gpvp->p[j]):(gpvp->p[j]>max1[j]?(max1[j]=gpvp->p[j]):i=i);
        }
    }
    gpvp=(Greedy_Projection_Vert_Prop*)(c2->vertices[0]->prop);
    double min2[2]={gpvp->p[0],gpvp->p[1] },
        max2[2]={gpvp->p[0],gpvp->p[1]};

    for(int i=1;i<c2->vertices_size;i++)
    {
        gpvp=(Greedy_Projection_Vert_Prop*)(c2->vertices[i]->prop);

        for(int j=0;j<2;j++)
        {
            gpvp->p[j]<min2[j]?(min2[j]=gpvp->p[j]):(gpvp->p[j]>max2[j]?(max2[j]=gpvp->p[j]):i=i);

        }
    }
    if(min2[0]>max1[0]+DEFAULT_THRESHOLD||min2[1]>max1[1]+DEFAULT_THRESHOLD)
    {
        return 0;
    }
    if(max2[0]<min1[0]+DEFAULT_THRESHOLD||max2[1]<min1[1]+DEFAULT_THRESHOLD)
    {
        return 0;
    }
    return 1;
}

static inline void touying_one_vertex(double*matd,template_v* v1)
{
    Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v1->prop);
    gpvp->is_init=1;
    gpvp->p[0]=(matd[0*3+0]*v1->point[0])+(matd[0*3+1]*v1->point[1])+(matd[0*3+2]*v1->point[2]);
    gpvp->p[1]=(matd[1*3+0]*v1->point[0])+(matd[1*3+1]*v1->point[1])+(matd[1*3+2]*v1->point[2]);
    gpvp->p[2]=(matd[2*3+0]*v1->point[0])+(matd[2*3+1]*v1->point[1])+(matd[2*3+2]*v1->point[2]);
    //printf("%lf %lf %lf\n",gpvp->p[0],gpvp->p[1],gpvp->p[2]);

}
static int delete_cell_is_manifold(Mesh * m,template_c* c )
{
    //return 0;
    Int_RB_Tree * tree_v=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree_v);

    //for(Node*nit=c->halffaces;nit!=NULL;nit=(Node*)(nit->Next))
    for(auto chfit=m->chf_begin(m,*c);chfit!=m->chf_end(m,*c);chfit++)
    {
        template_hf*hf=quote(chfit);
        if(Mesh_s_opposite_halfface(hf)->cell==NULL)
        {
            for(int i=0; i<hf->vertices_size;i++)
            {
                tree_v->insert(tree_v,hf->vertices[i]->id,hf->vertices[i]);
            }
        }
    }

    int re=1;
    for(int i=0;i<c->vertices_size;i++)
    {
        if(m->vertex_is_boundary(m,*(c->vertices[i])) )
        {
            if(tree_v->find(tree_v, c->vertices[i]->id)==NULL)
            {
                re=0;
                break;
            }
        }
    }

    int_rb_tree_free(tree_v);

    return re;
}
static double** my_get_shunxu_vertices_from_cell(template_c*c)
{
    double**re=(double**)malloc(sizeof(double*)*c->vertices_size);
    Node*temp_node= Mesh_adjust_halffaces(c);
    int i=0;
    for(Node* nit=temp_node;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_hf* hf=(template_hf*)(nit->value);
        Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(hf->vertices[1]->prop);
        re[i]=gpvp->p;
        i++;
    }

    free_node(temp_node);
    return re;
}
static int is_two_cells_intersecting2d(Mesh* m,template_c*c1,template_c*c2)
{
    for(auto chfit=m->chf_begin(m,*c1);chfit!=m->chf_end(m,*c1);chfit++)
    {

        if(Mesh_s_opposite_halfface(quote(chfit))->cell==c2)
        {
            template_f*f=quote(chfit)->face;

            // if(c1->vertices[0]->id==12756&&c1->vertices[1]->id==12757&&c1->vertices[2]->id==14935&&
            //     c2->vertices[0]->id==14935&&c2->vertices[1]->id==12756&&c2->vertices[2]->id==56183)
            // {
            //     dprintf("f:%d\n",f->id);
            //     dprintf("%d %d\n",f->halffaces[0]->vertices[0]->id,f->halffaces[1]->vertices[0]->id);
            // }
            //return 0;
            //return (f->halffaces[0]->vertices[0]==f->halffaces[1]->vertices[0]);
            return ((f->halffaces[0]->vertices[0]==f->halffaces[1]->vertices[0])?2:0);
        }
    }

    if(!is_two_cell_maybe_inter2d(c1,c2))
    {
        return 0;
    }

    int tuihua=0;
    double** poly=my_get_shunxu_vertices_from_cell(c1);
    // double **poly=(double**)malloc(sizeof(double*)*c1->vertices_size);
    // for(int i=0;i<c1->vertices_size;i++)
    // {
    //     Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(c1->vertices[i]->prop);
    //     poly[i]=gpvp->p;
    // }
    double** poly1=my_get_shunxu_vertices_from_cell(c2);
    // double **poly1=(double**)malloc(sizeof(double*)*c2->vertices_size);
    // for(int i=0;i<c2->vertices_size;i++)
    // {
    //     Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(c2->vertices[i]->prop);
    //     poly1[i]=gpvp->p;
    // }

    int len2=0;
    double**new_poly=(double**)malloc(sizeof(double*)*(c1->vertices_size+c2->vertices_size));
    for(int i=0;i<c1->vertices_size+c2->vertices_size;i++)
    {
        new_poly[i]=(double*)malloc(sizeof(double)*3);
    }
    get_intersections_from_two_polygons(poly,c1->vertices_size ,
        poly1 ,c2->vertices_size,new_poly,&len2);

    // double**new_poly= get_intersection_polygon_and_polygon(poly,c1->vertices_size,
    //     poly1,c2->vertices_size,&len2,&tuihua);

// void get_intersections_from_two_polygons(double**poly1,int len1,
//     double**poly2 ,int len2,double**result,int* len)

    free(poly);free(poly1);

    for(int i=0;i<c1->vertices_size+c2->vertices_size;i++)
    {
        free(new_poly[i]);
    }
    free(new_poly);


    return (len2!=0);
}

static int is_two_poly_intersecting2d(Mesh* m,template_c*c1,template_v* v1,template_v*v2,template_v*v3)
{
    for(auto chfit=m->chf_begin(m,*c1);chfit!=m->chf_end(m,*c1);chfit++)
    {
        if((quote(chfit)->vertices[1]==v1&& quote(chfit)->vertices[0]==v2)||
            (quote(chfit)->vertices[1]==v2&& quote(chfit)->vertices[0]==v3)||
            (quote(chfit)->vertices[1]==v3&& quote(chfit)->vertices[0]==v1))
        {
            return 0;
        }
    }

    int tuihua=0;
    double** poly=my_get_shunxu_vertices_from_cell(c1);
    // double **poly=(double**)malloc(sizeof(double*)*c1->vertices_size);
    // for(int i=0;i<c1->vertices_size;i++)
    // {
    //     Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(c1->vertices[i]->prop);
    //     poly[i]=gpvp->p;
    // }
    double*poly1[3]={((Greedy_Projection_Vert_Prop*)(v1->prop))->p,
        ((Greedy_Projection_Vert_Prop*)(v2->prop))->p,((Greedy_Projection_Vert_Prop*)(v3->prop))->p};


    int len2=0;
    // double**new_poly= get_intersection_polygon_and_polygon(poly,c1->vertices_size,poly1,3,&len2,&tuihua);
    double**new_poly=(double**)malloc(sizeof(double*)*(c1->vertices_size+3));
    for(int i=0;i<c1->vertices_size+3;i++)
    {
        new_poly[i]=(double*)malloc(sizeof(double)*3);
        memset(new_poly[i],0,sizeof(double)*3);
        // free(new_poly[i]);
    }

    get_intersections_from_two_polygons(poly,c1->vertices_size ,
        poly1 ,3,new_poly,&len2);

    free(poly);

    for(int i=0;i<c1->vertices_size+3;i++)
    {
        free(new_poly[i]);
    }
    free(new_poly);

    return 1;
}

static int is_two_poly_cell_intersecting3d(Mesh* m,template_c*c1,template_c*c2)
{
    for(auto chfit=m->chf_begin(m,*c1);chfit!=m->chf_end(m,*c1);chfit++)
    {
        if(Mesh_s_opposite_halfface(quote(chfit))->cell==c2)
        {
            return 0;
        }
    }
    if(is_two_cell_maybe_inter(c1,c2)==0)
    {
        return 0;
    }
    //double** poly1=my_get_shunxu_vertices_from_cell(c1);
    double **poly1=(double**)malloc(sizeof(double*)*c1->vertices_size);
    for(int i=0;i<c1->vertices_size;i++)
    {
        poly1[i]=c1->vertices[i]->point;
    }
    double **poly2=(double**)malloc(sizeof(double*)*c2->vertices_size);
    for(int i=0;i<c2->vertices_size;i++)
    {
        poly2[i]=c2->vertices[i]->point;
    }

    int len=0;
    double**mark_ps=(double**)malloc(sizeof(double*)*(c1->vertices_size+c2->vertices_size ));
    for(int i=0;i<c1->vertices_size+c2->vertices_size;i++)
    {
        mark_ps[i]=(double*)malloc(sizeof(double)*3);
        memset(mark_ps[i],0,sizeof(double)*3);
    }
    get_intersections_from_two_polygons(poly1,c1->vertices_size,
        poly2 ,c2->vertices_size,mark_ps ,&len);

    // double ** mark_ps=get_intersection_points_from_two_polygons(
    //     poly1,c1->vertices_size,poly2 ,c2->vertices_size,&jinhua);

    free(poly1);free(poly2);
    for(int i=0;i<c1->vertices_size+c2->vertices_size;i++)
    {
        free(mark_ps[i]);
    }
    free(mark_ps);
    return (len!=0);
}

static inline int is_triangle_clockwise(template_v*v1,template_v*v2,template_v*v3)
{
    Greedy_Projection_Vert_Prop* gpvp1=(Greedy_Projection_Vert_Prop*)(v1->prop);
    Greedy_Projection_Vert_Prop* gpvp2=(Greedy_Projection_Vert_Prop*)(v2->prop);
    Greedy_Projection_Vert_Prop* gpvp3=(Greedy_Projection_Vert_Prop*)(v3->prop);

    double p1[2]={gpvp1->p[0]-gpvp3->p[0],gpvp1->p[1]-gpvp3->p[1]};
    double p2[2]={gpvp2->p[0]-gpvp3->p[0],gpvp2->p[1]-gpvp3->p[1]};

    if(p1[0]*p2[1]-p1[1]*p2[0]>0)
    {
        return 1;
    }
    return 0;
}

static inline int is_cell_clockwise(template_c*c)
{
    template_hf* hf=(template_hf*)(c->halffaces[0]);

    template_v* v=NULL;
    for(int i=0;i<c->vertices_size;i++)
    {
        if(c->vertices[i]!=hf->vertices[0]&&c->vertices[i]!=hf->vertices[1])
        {
            v=c->vertices[i];
            break;
        }
    }
    if(v==NULL)
    {return 0;}
    return is_triangle_clockwise(hf->vertices[0],hf->vertices[1],v);
}
static Int_RB_Tree* filter_clockwise_cells(Int_RB_Tree*tree_c)
{
    Node* temp=NULL;
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c*c=(template_c*)(it.second);
        if(!is_cell_clockwise(c))
        {
            temp=node_overlying(temp,c);
        }
    }
    for(Node* nit=temp;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_c* c=(template_c*)(nit->value);
        tree_c->erase(tree_c,c->id);
    }


    free_node(temp);

    return tree_c;
}


static Int_RB_Tree* get_all_cells_of_vertices(Node* node_vs)
{
    Int_RB_Tree * re=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(re);
    for(Node* nit=node_vs;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v*v=(template_v*)(nit->value);
        // for(Node* nit1=v->cells;nit1!=NULL;nit1=(Node*)(nit1->Next))
        for(int i=0;i<v->cellsn.size;i++)
        {
            template_c* c=(template_c*)(v->cellsn.extra_data[i]);
            // template_c*c = (template_c*)(nit1->value);
            re->insert(re,c->id,c);
        }
    }

    return re;
}

static Int_RB_Tree* get_all_vertices_of_cells(Int_RB_Tree* tree_c)
{
    Int_RB_Tree * re=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(re);
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c*c=(template_c*)(it.second);
        for(int i=0;i<c->vertices_size;i++)
        {
            re->insert(re,c->vertices[i]->id,c->vertices[i]);
        }
    }

    return re;
}

static inline double* compute_triangle_normal(double* q1,double* q2,double* q3)
{
    double p1[3]={q2[0]-q1[0],q2[1]-q1[1],
        q2[2]-q1[2]};
    double p2[3]={q3[0]-q1[0],q3[1]-q1[1],
        q3[2]-q1[2]};
    double*re=out_product(p1,p2);
    normalize(re,3);
    return re;
}


static inline double compute_triangle_area(double* p1,double* p2,double* p3)
{
    double q1[3]={p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]};

    double q2[3]={p3[0]-p1[0],p3[1]-p1[1],p3[2]-p1[2]};

    double q[3]={q1[1]*q2[2]-q1[2]*q2[1],-q1[0]*q2[2]+q1[2]*q2[0],q1[0]*q2[1]-q1[1]*q2[0]};
    double re=SAFE_SQRT(q[0]*q[0]+q[1]*q[1]+q[2]*q[2])/2.0;

    return re;
}

static inline double compute_score_triangle(template_v* v1,template_v* v2,template_v* v3 )
{
    double normal[3]={0};
    Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v1->prop);
    normal[0]+=1.0*gpvp->normal[0]; normal[1]+=1.0*gpvp->normal[1];normal[2]+=1.0*gpvp->normal[2];
    gpvp=(Greedy_Projection_Vert_Prop*)(v2->prop);
    normal[0]+=gpvp->normal[0]; normal[1]+=gpvp->normal[1];normal[2]+=gpvp->normal[2];
    gpvp=(Greedy_Projection_Vert_Prop*)(v3->prop);
    normal[0]+=gpvp->normal[0]; normal[1]+=gpvp->normal[1];normal[2]+=gpvp->normal[2];
    normalize(normal,3);

    double* n1=compute_triangle_normal(v1->point,v2->point,v3->point );

    //double re=n1[0]*gpvp->normal[0]+n1[1]*gpvp->normal[1]+n1[2]*gpvp->normal[2]+1.0 ;
    double re=n1[0]*normal[0]+n1[1]*normal[1]+n1[2]*normal[2];

    SAFE_FREE(n1);

    return re;
}

static double* compute_cell_normal(template_c*c)
{
    if(c==NULL)
    {
        return NULL;
    }
    template_hf* hf=(template_hf*)(c->halffaces[0]);
    template_v* v=NULL;
    for(int i=0;i<c->vertices_size;i++)
    {
        if(c->vertices[i]!=hf->vertices[0]&&c->vertices[i]!=hf->vertices[1])
        {
            v=c->vertices[i];
            break;
        }
    }

    if(v==NULL)
    {return NULL;}

    double p1[3]={hf->vertices[0]->point[0]-v->point[0],hf->vertices[0]->point[1]-v->point[1],
        hf->vertices[0]->point[2]-v->point[2]};

    double p2[3]={hf->vertices[1]->point[0]-v->point[0],hf->vertices[1]->point[1]-v->point[1],
        hf->vertices[1]->point[2]-v->point[2]};

    double*re=out_product(p1,p2);
    normalize(re,3);
    return re;
}

static inline double compute_cell_area(template_c*c)
{
    if(c==NULL){return 0;}
    template_v* v=c->vertices[0];
    double re=0;
    //for(Node* nit=c->halffaces;nit!=NULL;nit=(Node*)(nit->Next))
    for(auto chfit=Mesh_chf_begin(NULL,*c);chfit!=Mesh_chf_end(NULL,*c);chfit++)
    {
        template_hf* hf=quote(chfit);
        if(hf->vertices[0]==v||hf->vertices[1]==v )
        {
            continue;
        }
        re+=compute_triangle_area(hf->vertices[0]->point,hf->vertices[1]->point,v->point);
    }

    return re;
}

static inline double compute_score_cell(template_c*c)
{
    double normal[3]={0};
    for(int i=0;i<c->vertices_size;i++)
    {
        Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(c->vertices[i]->prop);
        normal[0]+=1.0*gpvp->normal[0]; normal[1]+=1.0*gpvp->normal[1];normal[2]+=1.0*gpvp->normal[2];
    }
    normalize(normal,3);
    double*n1= compute_cell_normal(c);
    double re=n1[0]*normal[0]+n1[1]*normal[1]+n1[2]*normal[2];

    SAFE_FREE(n1);

    return re;
}


static  Node* get_qianzaidedian(Node* vs,double* n)
{
    Node* re=NULL;
    for(Node* nit=vs;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v*v=(template_v*)(nit->value);
        Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        if(inner_product(n,gpvp->normal,3)>0)
        {
            re=node_overlying(re,v);
        }
    }

    return re;
}

static Node* get_tree_c_active_boundary_hfs(Mesh* m,Int_RB_Tree* tree_c)
{
    Node* re=NULL;

    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c* c=(template_c*)(it.second);
        for(auto chfit=m->chf_begin(m,*c);chfit!=m->chf_end(m,*c);chfit++)
        {
            if(Mesh_s_opposite_halfface(quote(chfit))->cell==NULL)
            {
                Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(quote(chfit)->vertices[0]->prop);
                Greedy_Projection_Vert_Prop*gpvp1=(Greedy_Projection_Vert_Prop*)(quote(chfit)->vertices[1]->prop);
                if(gpvp->active==1&&gpvp1->active==1)
                {
                    re=node_overlying(re,quote(chfit));
                }
            }
        }
    }

    return re;
}

static int is_vertex_in_cell(template_v*v,template_c* c)
{
    //printf("begin dis vertex in cell\n");
    if(v==c->vertices[0]||v==c->vertices[1]||v==c->vertices[2])
    {
        return 0;
    }
    Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);

    template_hf* hf=(template_hf*)(c->halffaces[0]);

    Greedy_Projection_Vert_Prop*gpvp1=(Greedy_Projection_Vert_Prop*)(hf->vertices[0]->prop);
    Greedy_Projection_Vert_Prop*gpvp2=(Greedy_Projection_Vert_Prop*)(hf->vertices[1]->prop);
    double p1[2]={gpvp1->p[0]-gpvp->p[0],gpvp1->p[1]-gpvp->p[1]};
    double p2[2]={gpvp2->p[0]-gpvp->p[0],gpvp2->p[1]-gpvp->p[1] };
    double mark=p1[0]*p2[1]-p1[1]*p2[0];

    //for(Node* nit=c->halffaces;nit!=NULL;nit=(Node*)(nit->Next))

    for(auto chfit=Mesh_chf_begin(NULL,*c);chfit!=Mesh_chf_end(NULL,*c);chfit++)
    {
        hf=quote(chfit);

        gpvp1=(Greedy_Projection_Vert_Prop*)(hf->vertices[0]->prop);
        gpvp2=(Greedy_Projection_Vert_Prop*)(hf->vertices[1]->prop);
        p1[0]=gpvp1->p[0]-gpvp->p[0];p1[1]=gpvp1->p[1]-gpvp->p[1];
        p2[0]=gpvp2->p[0]-gpvp->p[0];p2[1]=gpvp2->p[1]-gpvp->p[1];
        if( (p1[0]*p2[1]-p1[1]*p2[0])*mark<=0)
        {
            return 0;
        }
    }
    //printf("en dis vertex in cell\n");
    return 1;
}
static Mesh* touying_chongjian_insert_vertex_to_cell(Mesh*m,template_v* v,template_c*c,Int_RB_Tree* tree_c)
{
    Node* temp=NULL;
    //printf("hrere1\n");
    //for(Node* nit=c->halffaces;nit!=NULL;nit=(Node*)(nit->Next))
    for(auto chfit=m->chf_begin(m,*c);chfit!=m->chf_end(m,*c);chfit++)
    {
        //printf("once %d %d ;",quote(chfit)->id,chfit.i);
        template_hf*hf=quote(chfit);
        if(compute_score_triangle(hf->vertices[0],hf->vertices[1],v )<=0)
        {
            free_node_value(temp);
            free_node(temp);
            return m;
        }
        template_v** vs=(template_v**)malloc(sizeof(template_v*)*3);
        vs[0]=hf->vertices[0];vs[1]=hf->vertices[1];vs[2]=v;
        temp=node_overlying(temp,vs);
    }

    tree_c->erase(tree_c,c->id);
    m->delete_cell(m,*c,true);
    for(Node* nit=temp;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v**vs=(template_v**)(nit->value);
        template_c*c1=m->create_cellv(m,vs,3);
        //if(c1!=NULL){tree_c->insert(tree_c,c1->id,c1);}
        tree_c->insert(tree_c,c1->id,c1);
        free(vs);
    }
    free_node(temp);

    return m;
}

static int is_triangle_inter_with_int_rb_tree(Mesh*m,Int_RB_Tree* tree_c,template_v* v1,template_v* v2,template_v*v3)
{
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c*c=(template_c*)(it.second);
        if(is_two_poly_intersecting2d(m,c,v1,v2,v3))
        {
            return 1;
        }
    }
    return 0;
}

static template_c* touying_chongjian_create_new_triangle2mesh(Mesh*m,Int_RB_Tree* tree_c,
    Node* node_vs)
{
    //printf("begin touying_chongjian_create_new_triangle2mesh");
    template_c* re=NULL;
    for(Node* nit=node_vs;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v*v1=(template_v*)(nit->value);
        if(nit->Next==NULL||v1->cellsn.size!=0){continue;}
        for(Node* nit1=(Node*)(nit->Next);nit1!=NULL;nit1=(Node*)(nit1->Next) )
        {
            template_v*v2=(template_v*)(nit1->value);

            if(nit1->Next==NULL||v2->cellsn.size!=0){continue;}

            for(Node* nit2=(Node*)(nit1->Next);nit2!=NULL;nit2=(Node*)(nit2->Next))
            {
                template_v*v3=(template_v*)(nit2->value);
                if(v3->cellsn.size!=0){continue;}
                if(!is_triangle_clockwise(v1,v2,v3))
                {
                    template_v* temp=v1;v1=v2;v2=temp;
                }
                if(compute_score_triangle(v1,v2,v3 )>0&&
                    !is_triangle_inter_with_int_rb_tree(m,tree_c,v1,v2,v3))
                {
                    template_v*vs[3]={v1,v2,v3};
                    re=m->create_cellv(m,vs,3);
                    tree_c->insert(tree_c,re->id,re);
                    return re;
                }
            }
        }

    }

    return re;
}


static template_c* touying_chongjian_insert_vertex(Mesh* m,template_v* v,Int_RB_Tree* tree_c,
    Node* node_vs,int is_insert_in)

{
    //printf("in touying_chongjian_insert_vertex:%d\n",tree_c->size);
    if(v->cellsn.size!=0){return NULL;}
    int i=0;
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        //printf("once %d ;",i);
        template_c* c=(template_c*)(it.second);
        if(is_vertex_in_cell(v,c))
        {
            if(is_insert_in==1)
            {
                 touying_chongjian_insert_vertex_to_cell(m,v,c,tree_c);
            }
            return NULL;
        }
        i++;
    }

    //printf("jrere1\n");
    Node*hfs= get_tree_c_active_boundary_hfs(m,tree_c);
    template_c*c1=NULL;
    for(Node* nit=hfs;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_hf*hf=(template_hf*)(nit->value);
        if(compute_score_triangle(hf->vertices[1],hf->vertices[0],v )>0&&
            is_triangle_clockwise(hf->vertices[1],hf->vertices[0],v)&&
            !is_triangle_inter_with_int_rb_tree(m,tree_c,hf->vertices[1],hf->vertices[0],v))
        {
            template_v*vs[3]={hf->vertices[1],hf->vertices[0],v};
            c1=m->create_cellv(m,vs,3);
            tree_c->insert(tree_c,c1->id,c1);

            break;
        }
    }
     //printf("jrere2\n");
    if(c1==NULL&&hfs==NULL)
    {
        c1=touying_chongjian_create_new_triangle2mesh(m,tree_c,node_vs);
    }
    free_node(hfs);

    return c1;
}
static int my_is_face_legal(Mesh*m,template_v*v1,template_v*v2)
{
    template_v*vs[2]={v1,v2};
    template_f*f=m->get_facev(m,vs,2);

    if(f==NULL){return 1;}

    if(m->face_is_boundary(m,*f) )
    {
        return 1;
    }


    return 0;
}
static void tiaozheng_tree_hfs(Mesh* m,Int_RB_Tree* tree_hfs,template_c* c)
{
    for(auto chfit=m->chf_begin(m,*c);chfit!=m->chf_end(m,*c);chfit++)
    {
        if(Mesh_s_opposite_halfface(quote(chfit) )->cell==NULL)
        {
            tree_hfs->insert(tree_hfs,quote(chfit)->id,quote(chfit) );
        }
        else
        {
            tree_hfs->erase(tree_hfs ,Mesh_s_opposite_halfface(quote(chfit))->id);
        }
    }

}
static int touying_chongjian_create_quadrilateral(Mesh* m,Int_RB_Tree* tree_hfs,Int_RB_Tree*tree_c)
{

    //dprintf("begin once quadri\n");
    for(auto it=tree_hfs->begin(tree_hfs);it.it!=NULL;it++)
    {
        template_hf*hf1=(template_hf*)(it.second);
        for(auto it1=it;it1.it!=NULL;it1++)
        {
            template_hf*hf2=(template_hf*)(it1.second);
            if(hf1==hf2){continue;}
            template_v* vs[2][3]={0};
            int temp_i=0;
            if(hf1->vertices[1]==hf2->vertices[0])
            {
                temp_i=1;vs[0][0]=hf1->vertices[1];vs[0][1]=hf1->vertices[0];vs[0][2]=hf2->vertices[1];
                if(!my_is_face_legal(m,hf1->vertices[0],hf2->vertices[1]))
                {
                    temp_i=0;
                }
            }
            else if(hf1->vertices[0]==hf2->vertices[1])
            {
                temp_i=1;vs[0][0]=hf1->vertices[1];vs[0][1]=hf1->vertices[0];vs[0][2]=hf2->vertices[0];
                if(!my_is_face_legal(m,hf1->vertices[1],hf2->vertices[0]))
                {
                    temp_i=0;
                }
            }
            else
            {
                temp_i=2;
                vs[0][0]=hf1->vertices[1];vs[0][1]=hf1->vertices[0];vs[0][2]=hf2->vertices[1];
                vs[1][0]=hf2->vertices[1];vs[1][1]=hf2->vertices[0];vs[1][2]=hf1->vertices[1];

                if(!my_is_face_legal(m,hf1->vertices[0],hf2->vertices[1])
                    ||!my_is_face_legal(m,hf1->vertices[1],hf2->vertices[0]))
                {
                    temp_i=0;
                }
            }
            //dprintf("hf1:%d %d hf2:%d %d temp_i:%d\n",hf1->vertices[0]->id,hf1->vertices[1]->id,hf2->vertices[0]->id,
            //    hf2->vertices[1]->id,temp_i);

            if(temp_i==1)
            {
                if(is_triangle_clockwise(vs[0][0],vs[0][1],vs[0][2])&&
                    compute_score_triangle(vs[0][0],vs[0][1],vs[0][2])>0&&
                    !is_triangle_inter_with_int_rb_tree(m,tree_c,vs[0][0],vs[0][1],vs[0][2])
                    )
                {
                    template_c*c=m->create_cellv(m,vs[0],3);
                    tiaozheng_tree_hfs(m,tree_hfs, c);
                    tree_c->insert(tree_c,c->id,c);
                    return 1;
                }
            }
            else if(temp_i==2)
            {
                if(is_triangle_clockwise(vs[0][0],vs[0][1],vs[0][2])&&
                    is_triangle_clockwise(vs[1][0],vs[1][1],vs[1][2]))
                {
                    if(!is_triangle_inter_with_int_rb_tree(m,tree_c,vs[0][0],vs[0][1],vs[0][2])&&
                        !is_triangle_inter_with_int_rb_tree(m,tree_c,vs[1][0],vs[1][1],vs[1][2]))
                    {
                        if(compute_score_triangle(vs[0][0],vs[0][1],vs[0][2])<=0||
                            compute_score_triangle(vs[1][0],vs[1][1],vs[1][2])<=0)
                        {
                            vs[0][0]=hf1->vertices[1];vs[0][1]=hf1->vertices[0];vs[0][2]=hf2->vertices[0];
                            vs[1][0]=hf2->vertices[1];vs[1][1]=hf2->vertices[0];vs[1][2]=hf1->vertices[0];
                        }

                        if(compute_score_triangle(vs[0][0],vs[0][1],vs[0][2])>0&&
                            compute_score_triangle(vs[1][0],vs[1][1],vs[1][2])>0&&
                            my_is_face_legal(m,vs[0][2],vs[1][2]))
                        {
                            template_c*c=m->create_cellv(m,vs[0],3);
                            tiaozheng_tree_hfs(m,tree_hfs, c);

                            tree_c->insert(tree_c,c->id,c);
                            c=m->create_cellv(m,vs[1],3);
                            tiaozheng_tree_hfs(m,tree_hfs, c);

                            tree_c->insert(tree_c,c->id,c);

                            return 1;
                        }
                    }
                }
            }

        }
    }
    return 0;
}

static Mesh* touying_chongjian_create_connected_halffaces(Mesh* m,Int_RB_Tree*tree_c )
{
    Node*hfs= get_tree_c_active_boundary_hfs(m,tree_c);

    Int_RB_Tree * tree_hfs=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree_hfs);
    for(Node* nit=hfs;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_hf*hf=(template_hf*)(nit->value);
        //dprintf("hf :%d %d\n",hf->vertices[0]->id,hf->vertices[1]->id);

        tree_hfs->insert(tree_hfs,hf->id,hf);
    }

    int flag= touying_chongjian_create_quadrilateral(m,tree_hfs,tree_c);

    while(flag)
    {
        flag= touying_chongjian_create_quadrilateral(m,tree_hfs,tree_c);
    }



    int_rb_tree_free(tree_hfs);
    free_node(hfs);
    return m;
}

//***************
//大于0是在外心圆内部
//小于0是在外心圆外部
//****************
static inline double is_concyclic_four_vertices(double* p1,double* p2,double* p3,double* p4)
{
    LB_Matrix* mat=(LB_Matrix*)malloc(sizeof(LB_Matrix));
    lb_matrix_init_double(mat,4,4);
    double* data=(double*)(mat->data);
    data[0*4+0]=p1[0]; data[0*4+1]=p1[1];data[0*4+2]=p1[0]*p1[0]+p1[1]*p1[1];data[0*4+3]=1.0;
    data[1*4+0]=p2[0]; data[1*4+1]=p2[1];data[1*4+2]=p2[0]*p2[0]+p2[1]*p2[1];data[1*4+3]=1.0;
    data[2*4+0]=p3[0]; data[2*4+1]=p3[1];data[2*4+2]=p3[0]*p3[0]+p3[1]*p3[1];data[2*4+3]=1.0;
    data[3*4+0]=p4[0]; data[3*4+1]=p4[1];data[3*4+2]=p4[0]*p4[0]+p4[1]*p4[1];data[3*4+3]=1.0;
    double *re=((double*)(mat->det(mat)));
    double re1=*re;
    lb_matrix_free(mat);
    free(re);
    return re1;
}

static inline double compute_weight_score_2triangles(template_v* v1,template_v* v2,template_v* v3,
    template_v* v4,template_v* v5,template_v* v6 )
{
    double score1=compute_score_triangle(v1,v2,v3);
    double weight1=compute_triangle_area(v1->point,v2->point,v3->point);
    double score2=compute_score_triangle(v4,v5,v6);
    double weight2=compute_triangle_area(v4->point,v5->point,v6->point);

    return (score1*weight1+score2*weight2)/(weight1+weight2);
}

static inline  template_f* swap_edge(Mesh* m, template_f* f,template_v** vs,Int_RB_Tree* tree_c)
{
    //printf("begin swap\n");

    template_c* c1=f->halffaces[0]->cell;
    template_c* c2=f->halffaces[1]->cell;
    if(c1==NULL||c2==NULL){return NULL;}



    tree_c->erase(tree_c,c1->id);
    tree_c->erase(tree_c,c2->id);

    m->delete_cell(m,*c1,true);
    m->delete_cell(m,*c2,true);

    template_v* vs1[3]={vs[3],vs[0],vs[1]};

    c1=m->create_cellv(m,vs1,3);
    tree_c->insert(tree_c,c1->id,c1);

    vs1[0]=vs[1];vs1[1]=vs[2];vs1[2]=vs[3];
    c2=m->create_cellv(m,vs1,3);
    tree_c->insert(tree_c,c2->id,c2);

    template_v* vs2[2]={vs[1],vs[3]};
    template_f* re= m->get_facev(m,vs2,2);

    return re;
}

static template_f* touying_tiaozheng_edge(Mesh* m, template_f* f,Int_RB_Tree* tree_c)
{
    template_v* vs[4]={0};

    template_c* c1=f->halffaces[0]->cell;
    template_c* c2=f->halffaces[1]->cell;
    if(c1==NULL||c2==NULL||c1->vertices_size>3||c2->vertices_size>3){return NULL;}
    //if(c1==NULL||c2==NULL){return NULL;}
    vs[0]=f->halffaces[0]->vertices[0]; vs[2]= f->halffaces[0]->vertices[1];
    for(int i=0;i<3;i++)
    {
        template_v* v=c1->vertices[i];
        if(v!=vs[0]&&v!=vs[2])
        {
            vs[3]=v;
            break;
        }
    }
    for(int i=0;i<3;i++)
    {
        template_v* v=c2->vertices[i];
        if(v!=vs[0]&&v!=vs[2])
        {
            vs[1]=v;
            break;
        }
    }
    template_v*vs1[2]={vs[1],vs[3]};
    if(m->get_facev(m,vs1,2)!=NULL){return  NULL;}
    Greedy_Projection_Vert_Prop*gpvps[4]={0};
    gpvps[0]=(Greedy_Projection_Vert_Prop*)(vs[0]->prop);
    gpvps[1]=(Greedy_Projection_Vert_Prop*)(vs[1]->prop);
    gpvps[2]=(Greedy_Projection_Vert_Prop*)(vs[2]->prop);
    gpvps[3]=(Greedy_Projection_Vert_Prop*)(vs[3]->prop);

    if(is_concyclic_four_vertices(gpvps[0]->p,gpvps[1]->p,gpvps[2]->p,gpvps[3]->p)<=0)
    {
        return NULL;
    }

    double score0=compute_weight_score_2triangles(vs[1],vs[2],vs[0],
        vs[3],vs[0],vs[2]);

    double score1=compute_score_triangle(vs[0],vs[1],vs[3]);
    double weight1=compute_triangle_area(vs[0]->point,vs[1]->point,vs[3]->point);
    double score2=compute_score_triangle(vs[2],vs[3],vs[1]);
    double weight2=compute_triangle_area(vs[2]->point,vs[3]->point,vs[1]->point);

    // double score2=compute_weight_score_2triangles(vs[0],vs[1],vs[3],
    //     vs[2],vs[3],vs[1]);

    //printf("score1:%lf score2:%lf\n",score1,score2);

    if(score1>0&&score2>0&& ( (score1*weight1+score2*weight2)/(weight1+weight2))>=1.0*score0)
    {
        return swap_edge(m,f,vs,tree_c);
    }
    return NULL;
}


static inline int my_panduan_f_legal(template_f*f,int *mark)
{
    return (f->halffaces[0]->cell->prop==mark&&f->halffaces[1]->cell->prop==mark);
}
static void tiaozheng_tree_cellsn(Mesh* m,Int_RB_Tree* tree_c)
{
    Int_RB_Tree * tree_f=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree_f);
    int mark=0;
    void** temp_pro=(void**)malloc(sizeof(void*)*tree_c->size);
    int i=0;
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c* c=(template_c*)(it.second);
        temp_pro[i]=c->prop;
        c->prop=&mark;
        i++;
    }

    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c* c=(template_c*)(it.second);
        for(auto chfit=m->chf_begin(m,*c);chfit!=m->chf_end(m,*c);chfit++)
        {
            if(!m->face_is_boundary(m,*(quote(chfit)->face)))
            {
                if(my_panduan_f_legal(quote(chfit)->face,&mark))
                {
                    tree_f->insert(tree_f,quote(chfit)->face->id,quote(chfit)->face);
                }
            }
        }
    }

    i=0;
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c* c=(template_c*)(it.second);
        c->prop=temp_pro[i];
        i++;
    }
    free(temp_pro);
    int flag=1;
    int num=tree_f->size;
    while(flag)
    {
        Int_RB_Tree * tree_f1=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
        int_rb_tree_init(tree_f1);
        flag=0;
        for(auto it=tree_f->begin(tree_f);it.it!=NULL;it++)
        {
            template_f*f=m->get_facep(m,it.first);
            if(f==NULL||m->face_is_boundary(m,*f))
            {
                continue;
            }
            template_f*f1= touying_tiaozheng_edge(m, f,tree_c);
            if(f1!=NULL)
            {
                tree_f1->insert(tree_f1,f1->id,f1);flag=1;
                flag=1;
            }
            else
            {
                f=m->get_facep(m,it.first);
                if(f!=NULL)
                {
                    tree_f1->insert(tree_f1,f->id,f);
                }
            }
        }
        int_rb_tree_free(tree_f);
        tree_f=tree_f1;
        if(num<=1){break;}
        //printf("tiaozheng\n");
        num--;
    }

    int_rb_tree_free(tree_f);
}

static Mesh* delete_intersecting_cells2d(Mesh*m ,Int_RB_Tree*tree_c)
{
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c*c=(template_c*)(it.second);
        //dprintf("c :%d %d %d\n",c->vertices[0]->id,c->vertices[1]->id,c->vertices[2]->id);

    }
    //if(!global_debug){return m;}
    int len=tree_c->size;
    template_c**cells=(template_c**)malloc(sizeof(template_c*)*tree_c->size);
    memset(cells,0,sizeof(template_c*)*tree_c->size);
    int** marks=(int**)malloc(sizeof(int*)*tree_c->size);
    int i=0;
    for(auto it =tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c* c=(template_c*)(it.second);
        cells[i]=c;
        marks[i]=(int*)malloc(sizeof(int)*tree_c->size);
        memset( marks[i],0,sizeof(int)*tree_c->size);
        i++;
    }

    int flag=1;
    while(flag)
    {
        flag=0;
        for(i=0;i<len;i++)
        {
            if(cells[i]==NULL){continue;}
            for(int j=i+1;j<len;j++)
            {
                if(cells[j]==NULL||marks[i][j]==-1){continue;}
                if(marks[i][j]==0)
                {
                    marks[i][j]=is_two_cells_intersecting2d(m,cells[i],cells[j]);
                    if(marks[i][j]==0)
                    {
                        marks[i][j]=-1;
                    }
                }
                int index=-1;
                if(marks[i][j]>=1 )
                {
                    if(delete_cell_is_manifold(m,cells[i]))
                    {
                        index=i;;
                        if(delete_cell_is_manifold(m,cells[j]))
                        {
                            if(compute_score_cell(cells[i])>compute_score_cell(cells[j]))
                            {index=j;}
                        }
                    }
                    else
                    {
                        if(delete_cell_is_manifold(m,cells[j]))
                        {index=j;}
                        else if(marks[i][j]==2){
                            index=i;
                            if(compute_score_cell(cells[i])>compute_score_cell(cells[j]))
                            {index=j;}
                        }
                    }
                    marks[i][j]=-1;
                }

                if(index!=-1)
                {
                    tree_c->erase(tree_c,cells[index]->id);
                    m->delete_cell(m,*(cells[index]),true);
                    cells[index]=NULL;
                    flag=1;
                    break;
                }
            }
            if(flag){break;}
        }
    }
    free(cells);
    for(i=0;i<len;i++)
    {
        free(marks[i]);
    }
    free(marks);
    return m;
}
Mesh* projective_reconstruct_fix(Mesh* m,Node* node_vs,double* n)
{
    Int_RB_Tree* tree_c=get_all_cells_of_vertices(node_vs);
    //对tree_c删除3d交叉面片//

    // double theta = acos(n[2]);
    // double  n1[3]={n[1],-n[0],0};
    // double*matd = get_matrix_from_rotation(theta,n1);
    double n1[3]={0,0,1};
    double *matd=get_one_dir_trans2dir_rota_matrix(n,n1);
    Node* node_vs1 = get_qianzaidedian(node_vs,n);
    for(Node* nit=node_vs1;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v* v=(template_v*)(nit->value);
        //dprintf("node_vs1:%d\n",v->id);
        Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        gpvp->is_init=1;gpvp->active=1;
        touying_one_vertex(matd,v);
    }
    Int_RB_Tree* tree_c1=get_all_cells_of_vertices(node_vs1);

    Int_RB_Tree*tree_v= get_all_vertices_of_cells(tree_c1);

    for(auto it=tree_v->begin(tree_v);it.it!=NULL;it++)
    {
        template_v*v=(template_v*)(it.second);
        Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        if(gpvp->is_init==0)
        {
            touying_one_vertex(matd,v);
        }
    }

    filter_clockwise_cells(tree_c1);


    // 对于tree_c1删除2d交叉面片 //

    delete_intersecting_cells2d(m,tree_c1);

//
    for(Node* nit=node_vs1;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v* v=(template_v*)(nit->value);
        Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        gpvp->is_init=0;gpvp->active=0;
    }
    for(auto it=tree_v->begin(tree_v);it.it!=NULL;it++)
    {
        template_v*v=(template_v*)(it.second);
        Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        gpvp->is_init=0;
    }
    free(matd);
    int_rb_tree_free(tree_v);
    int_rb_tree_free(tree_c1);
    free_node(node_vs1);
    int_rb_tree_free(tree_c);

    return m;
}

Mesh* projective_reconstruct(Mesh* m,Node* node_vs,double* n,int is_insert_in)

{
    Int_RB_Tree* tree_c=get_all_cells_of_vertices(node_vs);
    //对tree_c删除3d交叉面片//

    // double theta = acos(n[2]);
    // double  n1[3]={n[1],-n[0],0};
    // double*matd = get_matrix_from_rotation(theta,n1);
    double n1[3]={0,0,1};
    double *matd=get_one_dir_trans2dir_rota_matrix(n,n1);
    Node* node_vs1 = get_qianzaidedian(node_vs,n);
    for(Node* nit=node_vs1;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v* v=(template_v*)(nit->value);
        //dprintf("node_vs1:%d\n",v->id);
        Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        gpvp->is_init=1;gpvp->active=1;
        touying_one_vertex(matd,v);
    }
    //dprintf("node_vs1 size:%d\n",node_size(node_vs1));

    //dprintf("here1\n");
    Int_RB_Tree* tree_c1=get_all_cells_of_vertices(node_vs1);

    Int_RB_Tree*tree_v= get_all_vertices_of_cells(tree_c1);

    for(auto it=tree_v->begin(tree_v);it.it!=NULL;it++)
    {
        template_v*v=(template_v*)(it.second);
        Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        if(gpvp->is_init==0)
        {
            touying_one_vertex(matd,v);
        }
    }
    //printf("begin  filter_clockwise_cells\n");
    filter_clockwise_cells(tree_c1);

    //对于tree_c1删除2d交叉面片//
    //delete_intersecting_cells2d(m,tree_c1);
    //printf("begin  touying_chongjian_insert_vertex :%d\n",node_size(node_vs1));
    for(Node* nit=node_vs1;nit!=NULL;nit=(Node*)(nit->Next))
    {

        touying_chongjian_insert_vertex(m,(template_v*)(nit->value),tree_c1,node_vs1,is_insert_in);

    }
    //dprintf("begin connected halffaces\n");
    //dprintf("tree_c1 size:%d\n",tree_c1->size);
    //printf("begin create connetex\n");
    //if(global_debug==0)
    {
        touying_chongjian_create_connected_halffaces(m,tree_c1 );
    }
    //printf("end create connetex\n");

    //if(global_debug==0)
    {
        tiaozheng_tree_cellsn(m,tree_c1);

    }

//
    //dprintf("node_vs1 size:%d\n",node_size(node_vs1));

    for(Node* nit=node_vs1;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v* v=(template_v*)(nit->value);
        Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        gpvp->is_init=0;gpvp->active=0;
    }
    for(auto it=tree_v->begin(tree_v);it.it!=NULL;it++)
    {
        template_v*v=(template_v*)(it.second);
        Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        gpvp->is_init=0;
    }
    free(matd);
    int_rb_tree_free(tree_v);
    int_rb_tree_free(tree_c1);
    free_node(node_vs1);
    int_rb_tree_free(tree_c);

    return m;
}
static template_c* test_create_hole_triangle(Mesh* m, template_v* v,int fix,Node* node_v)
{
    template_v* v1 = NULL, * v2 = NULL;
    for (auto vfit = m->vf_begin(m, *v); vfit != m->vf_end(m, *v); vfit++)
    {
        if (m->face_is_boundary(m, *vfit))
        {
            template_hf* hf = NULL;
            if (quote(vfit)->halffaces[0]->cell != NULL)
            {
                hf = quote(vfit)->halffaces[0];
            }
            else
            {
                hf = quote(vfit)->halffaces[1];
            }
            if (hf->cell == NULL) { printf("is  null\n");return NULL; }
            if (hf->vertices[0] == v) { v1 = hf->vertices[1]; }
            else { v2 = hf->vertices[0]; }
        }
    }
    //compute_weight_score_triangle
    //compute_score_triangle

    if (v1==NULL||v2==NULL)
    {
        return NULL;
    }

    template_v* vs1[2] = { v1,v2 };
    template_f* f = m->get_facev(m, vs1, 2);
    if (f != NULL && !m->face_is_boundary(m, *f))
    {
        return NULL;
    }
    double*n= compute_triangle_normal(v1->point,v->point,v2->point);
    for(Node* nit=node_v;(nit!=NULL&&fix!=1);nit=(Node*)(nit->Next))
    {
        template_v* vv=(template_v*)(nit->value);
        if(vv== v||vv==v1||vv==v2){continue;}
        if(vv->cellsn.size != 0 && m->vertex_is_boundary(m, *vv))
        {
            double* nn1=compute_triangle_normal(vv->point,v1->point,v->point);
            double* nn2=compute_triangle_normal(vv->point,v->point,v2->point);
            double* nn3=compute_triangle_normal(vv->point,v2->point,v1->point);

            if(inner_product(nn1,n,3)>0&&inner_product(nn2,n,3)>0&&
                inner_product(nn3,n,3)>0)
            {
                free(n);free(nn1);free(nn2);free(nn3);
                return NULL;
            }
            free(nn1);free(nn2);free(nn3);
        }
    }
    SAFE_FREE(n);
    if (fix==1||compute_score_triangle(v1, v, v2) > 0)
    {
        int num = 0;
        for (auto vfit = m->vf_begin(m, *v1); vfit != m->vf_end(m, *v1); vfit++)
        {
            if (m->face_is_boundary(m, *vfit))
            {
                num++;
            }
        }

        template_v* vs[3] = { v1,v,v2 };
        template_c*cc= m->create_cellv(m, vs, 3);
        return cc;
    }
    return NULL;
}

void projective_reconstruct_fix_neighborhood_hole(Mesh* m, Node* node_v, double* normal,int fix)
{
    for (Node* nit = node_v; nit != NULL; nit = (Node*)(nit->Next))
    {
        template_v* v = (template_v*)(nit->value);

        if (v->cellsn.size != 0 && m->vertex_is_boundary(m, *v)&&m->is_manifold_vertices(m, v) )
        {

            test_create_hole_triangle(m, v,fix,node_v);
            //free_n_cells_node(n_cells);
        }

    }

}
